## etl -- 2.2.1 An Explicit Tau-Leaping Method

etl <- function(tfinal, x0, V, transitions)
{
	Vn <- V
	Vn[Vn>=0] <- 0
	Vn <- abs(Vn)
	nc <- 10
	epsilon <- 0.03

	df <- data.frame()
	nextstep <- 0
	x <- x0
	t <- 0
	while(t < tfinal) {
		## step 1 -- identify all critical transitions
		if(nextstep <= 1) {
			L <- floor(apply(x0/Vn, 2, min))
			Jcr <- L < nc
			Jncr <- L >= nc
			nextstep <- 0
		}
		## step 2 -- bound tau by epsilon
		if(nextstep <= 2) {
			a <- transitions(x)
			asum <- sum(a)
			means <- apply(V[,Jcr]*a[Jcr], 1, sum)
			vars <- apply((V[,Jcr]^2)*a[Jcr], 1, sum)
			g <- 1  # TODO (see Cao): fortunately, in our case all is 1st order
			maxs <- sapply(epsilon*x/g, function(x) max(x,1))
			tau1 <- min(maxs/abs(means), maxs/vars)
			nextstep <- 0
		}
		## step 3 -- do 100 SSA if tau<10*(1/lambda), then go back
		if(nextstep <= 3) {
			if(tau1 < (1/asum)*10) {
				for(n in 1:100) {
					a <- transitions(x)
					asum <- sum(a)
					tau <- rexp(1, asum)
					i <- which(runif(1) < cumsum(a)/asum)[1]
					t <- t + tau
					x <- x + V[,i]
					df <- rbind(df, data.frame(time=t, U=x[1], C=x[2], J=x[3]))
				}
				nextstep <- 2
			}
			else
				nextstep <- 0
		}
		## step 4 -- generate a second candidate time leap
		if(nextstep <= 4) {
			lambdac <- sum(a[Jcr])
			tau2 <- rexp(1, lambdac)
			nextstep <- 0
		}
		## step 5 -- set tau, and approximate K number of transitions
		if(nextstep <= 5) {
			tau <- min(tau1, tau2)
			K <- rep(0, ncol(V))
			K[Jncr] <- rpois(1, 1/(a[Jncr]*tau))
			if(tau==tau2) {
				j <- which(runif(1) < cumsum(a[Jcr])/lambdac)[1]
				K[which(Jcr)[j]] <- 1
			}
			nextstep <- 0
		}
		## step 6 -- 
		if(nextstep <= 6) {
			x_ <- x + apply(K*t(V),2,sum)
			if(any(x_ < 0)) {
				tau1 <- tau1/2
				nextstep <- 3
			}
			else {
				t <- t + tau
				x <- x_
				nextstep <- 0
				df <- rbind(df, data.frame(time=t, U=x[1], C=x[2], J=x[3]))
			}
		}
	}
	df
}
